NetNews Offline 2/NetNews Offline Volume 2.iso/news/comp/sys/amiga/applications/560

Text File | 1996-08-05 | 28.8 KB | 801 lines

Path: news.restena.lu!usenet
From: [e-mail] (Manou BILLA)
Newsgroups: comp.sys.amiga.applications
Subject: Re: computer algebra prog on Amiga  ?
Date: 13 Jan 1996 02:53:36 GMT
Organization: Not organized
Message-ID: <[e-mail]>
References: <4d1k6m$[e-mail]>
Reply-To: [e-mail]
NNTP-Posting-Host: slip4.restena.lu
X-Newsreader: THOR 2.22 (Amiga;TCP/IP)


On 11-Jan-96 01:04:38, COLET VINCENT ([e-mail]) wrote to All () about
computer algebra prog on Amiga  ? (<4d1k6m$[e-mail]>):

Hi COLET VINCENT

>Can we found "computer algebra" handle programmes on Amiga ?
>(like Mathematica, Mapple, ... on PC/Mac ...)

Yes there's MapleV R3 for the AMIGA. I don't know the price but it seems to be
worth it! (I don't have it :-(()

But there's a program called 'gp amiga' (plus pari though not needed).
This program is PD. You can't do as much symbolic math as Mathematica,
Maple, Reduce ro MacSyma with it, but this program is unbeatable when doing
numerical maths!  with Pari symbolic math is possible, but you can't compare
it to one of the above programs.

It is very fast. According to the doc it does compute up to 5 - 100 times
faster than one of the programs mentioned above.

The current gp (Amiga version) is version 1.38.3 (The on I use) comes in
several compiled versions (depends on CPU used).
The program is a port from UN*X boxes (SPARC, HP, VAX, MAC, PC, ATARI ...).
It is SHELL driven but you can even plot graphs and output them to postscript
printers and files or as LaTeX files.
The included DOC (in DVI format) is about 150 DIN A4 pages with a lot of
examples and explanations of the PARI programming language (much like LISP).

The program is comparable to musimp/mumath on CP/M computers some years ago.

Here's the AMIGA.readme:

About PARI             (14 dec 1993)
----------

PARI is the name of a sophisticated and free math package. GP is a
calculator that offers all the features of PARI and some more. PARI
uses *infinite* precision rational numbers and *arbitrary* precision
floating point numbers.

You can use complex numbers, vectors, matrices, polynomials, rational
functions and taylor expansions. PARI also handles integers mod n,
finite fields, algebraic numbers and p-adic numbers. PARI includes
standard numerical methods and the GP calculator also includes
hi-resolution plotting.

PARI is written by four professional number theorists, C. Batut,
D. Bernardi, H. Cohen and M. Olivier. The latter two are Professors of
Mathematics.

-------

This amiga distribution contains the GP calculator compiled for
different processors, an emacs mode for running GP, partial
documentation, and all the amiga specific files I used to compile GP.
The amiga hi-resolution plotting functions are written by Jerry
Tunnell, who kindly let me use and distribute them.

If you want full source and documentation, you will have to get the
source distribution. It *should* be available where you found this
package, as file pari-1.38.3.gz or something similar (*Please* keep
this archive and the source archive together. If you like GP, you will
probably want the documentation too). If you can't find the source
anywhere else, you can try to ftp to megrez.ceremab.u-bordeaux.fr,
directory pub/pari/unix. This is the main PARI site.


Files
-----

amiga/                   Amiga specific files and sources.
        makefile.68000   Makefiles for different amiga versions.
        makefile.68020
        makefile.68881
        mpAmiga.s        Assembler file (gcc syntax) for the 68020 versions.
                         Converted from mp.s with the convert68k.el program in
                         the elisp directory.
        plotAmiga.c      Hi-resolution plotting functions, written by J.B.
Tunnell.
        version68k.diff  Source diffs to add an Amiga version string.
        versionport.diff

bin/
        gp.68000         GP binaries for different processors.
        gp.68020
        gp.68881

doc/                     This directory does not contain the complete
documentation,
        usersch3.tex     only one file that is needed by pari.el

elisp/
        convert68k.el    Elisp program to convert a sun3 style 68k assembler
file
                         (read mp.s) into something that amiga gcc can
understand.
        pari.el          An Emacs mode for the GP calculator. Desribed below
and
        pari.elc         in the file pari.txt
        pari.menu        Used by pari.el.
        pari.txt         A description of pari-mode

examples/
        EXPLAIN          Description of the examples.
        Makefile         Note that you cannot compile the C example
        Makesimple       without the libpari.a library.
        bench.gp
        clareg.gp
        lucas.gp
        mattrans.c
        rho.gp
        squfof.gp
        tutnf.gp
        tutnfout


Starting PARI
-------------

First, you need to install Markus Wild's ixemul.library, if you don't
have it already. Version 39.45 is the most recent non-buggy version I
know of (39.47 seems to be unreliable). This library is available on
Aminet (for example at ftp.luth.se) and is included in the gcc
distribution.

GP (file gp.68020 or whichever version you use) takes three command line
options. The most important is '-s STACKSIZE'. This sets the initial
size of the internal PARI stack (not to be confused with the task
stack). The default value is 4 MB which may be more RAM than you have
available. Try 'gp -s 1000000' or 'gp -s 100000' if GP refuses to
start. The other two flags are '-p PRIMELIMIT' and '-b BUFFERSIZE'.
Default values are 500000 and 30000 respectively.

Talking about the task stack, I don't know exactly how large it must
be. I use a stack of 100000 bytes, and that seems to be enough. To set
the task stack, use the command 'STACK 100000' command from the shell,
not the -s option to GP.

At the pari command prompt (default '?'), \q or CTRL-\ exits GP. You
can type '?' to get some on-line help. Note that running GP inside
emacs gives you better online help.

The GP command interface is quite straight forward if you are used to
MATLAB or similar systems. Note that with GP both vectors and matrices
are typed with with square brackets '[' ']', with comma ',' separating
elements on the same row and semicolon ';' separating rows. For
example, a 2-2 matrix is typed '[1,2 ; 3,4]'.

The emacs mode.
---------------
To use this on the amiga, you must make sure that you have mounted the
FIFO: device, and that the SHELL environment variable is set to some
unix-style shell. I use the shell distributed with gcc, a port of
pdksh (file name gcc/bin/sh). The shell distributed with GNUemacs
might work too, but I haven't tried it. The emacs mode is described in
the file elisp/pari.txt. If you don't wan't to edit the pari.el file,
you should assign PARI: to the directory where you have installed
PARI.

Known bugs
----------
GP does not respond to CTRL-C when run from the shell. However, if you
send the CTRL-C signal from another shell window (with the BREAK
command) or type CTRL-C in GP's emacs buffer, GP is interrupted.

A free() call occasionally failes when using the 68000 or 68020
versions of GP. I have not had this problem with the 68881 version.

I'm tempted to blame both these problems on the ixemul.library, but
I'm not sure what happens.

For those who are curious about the differences between the three
versions gp.68000, gp.68020 and gp.68881:

* The first two are compiled with gcc -msoftfloat instead of gcc
-m68881. If a 68881 processor is present, all three version makes use
of it. The performance difference between the gp.68020 and the 68881
version should be rather small on any machine that can run both. I
included the 68881 because it seemes more reliable.

* In the 68020 and 68881 versions, some low level functions are
written in 68020 assembler, while the 68000 version is written
entirely in C and is compiled with gcc -m68000 to make sure that it
contains only 68000 instructions.

-------
Enjoy GP!

Feel free to send me comments and questions (and even bug reports).
        Niels M÷ller
        StΣlldalsvΣgen 11
        122 43 Enskede
        SWEDEN
email:  [e-mail]

For questions and bug reports not specific to the amiga version, you
can also write to the authors:
        [e-mail]


--------------------------- end readme AMIGA: gp amiga/ pari ------


And here's an example of the calculation of 1000! (factorial of 1000)

Calculation time with display of the result on my A3000 (68030, 68882, MMU, 25
MHz, with 12 MB Fast, 2 MB CHIP) on a 1024x756x8 Workbench screen (TIGA gfx
board) was:

                3 secs!!!


Here's the result: (sorry for the bad formatting of the output, but that is
                    the fault of my text editor)

Work2:science/gp/gp.68881
            GP/PARI CALCULATOR Version 1.38.3
                      (68020 version)

Copyright 1989,1993 by C. Batut, D. Bernardi, H. Cohen and M. Olivier

Type ? for help

\precision      = 28
\serieslength   = 16
\format         = g0.28
\prompt         = ?
stacksize = 4000000, prime limit = 500000, buffersize = 30000
? ?

1: Standard monadic or dyadic OPERATORS
2: CONVERSIONS and similar elementary functions
3: TRANSCENDENTAL functions
4: NUMBER THEORETICAL functions
5: Functions related to ELLIPTIC CURVES
6: Functions related to general NUMBER FIELDS
7: POLYNOMIALS and power series
8: Vectors, matrices and LINEAR ALGEBRA
9: SUMS, products, integrals and similar functions
10: GRAPHIC functions
11: PROGRAMMING under GP

Further help: ?n (1<=n<=11), ?functionname, or ?\ (keyboard commands)

? 1000!
%1 =
402387260077093773543702433923003985719374864210714632543799910429938512398629
020592044208486969404800479988610197196058
631666872994808558901323829669944590997424504087073759918823627727188732519779
50595099527612087497546249704360141827809464649
629105639388743788648733711918104582578364784997701247663288983595573543251318
53239584630755574091142624174743493475534286465
766116677973966688202912073791438537195882498081268678383745597317461360853795
34524221586593201928090878297308431392844403281
231558611036976801357304216168747609675871348312025478589320767169132448426236
13141250878020800026168315102734182797770478463
586817016436502415369139828126481021309276124489635992870511496497541990934222
15668325720808213331861168115536158365469840467
089756029009505376164758477284218896796462449451607653534081989013854424879849
59953319101723355556602139450399736280750137837
615307127761926849034352625200015888535147331611702103968175921510907788019393
17811419454525722386554146106289218796022383897
147608850627686296714667469756291123408243920816015378088989396451826324367161
67621791689097799119037540312746222899880051954
444142820121873617459926429565817466283029555702990243241531816172104658320367
86906117260158783520751516284225540265170483304
226143974286933061690897968482590125458327168226458066526769958652682272807075
78139185817888965220816434834482599326604336766
017699961283186078838615027946595513115655203609398818061213855860030143569452
72242063446317974605946825731037900840244324384
656572450144028218852524709351906209290231364932734975655139587205596542287497
74011413346962715422845862377387538230483865688
976461927383814900140767310446640259899490222221765904339901886018566526485061
79970235619389701786004081188972991831102117122
984590164192106888438712185564612496079872290851929681937238864261483965738229
11231250241866493531439701374285319266498753372
189406942814341185201580141233448280150513996942901534830776445690990731524332
78288269864602789864321139083506217095002597389
863554277196742822248757586765752344220207573630569498825087968928162753848863
39690995982628095612145099487170124451646126037
902930912088908694202851064018215439945715680594187274899809425474217358240106
36774045957417851608292301353580818400969963725
242305608559037006242712434169090041536901059339838357779394109700277534720000
00000000000000000000000000000000000000000000000
000000000000000000000000000000000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000000
?

---------------------- end example 1000! ----------------------

Here's the readme text of the included examples:

Several examples of complete and non-trivial GP programs are given in this
directory (as well as the C program mattrans.c using the Pari library
described
in Chapter 4 of the users' manual). This file gives a brief description of
these programs. All these programs should be read into GP by the command
\r file.

1) bench.gp

This program computes the first 1000 terms of the Fibonacci sequence, the
product p of successive terms, and the lowest common multiple q. It outputs
the ratio log(p)/log(q) every 10 terms (this ratio tends to pi^2/6 as k
tends to infinity). The name bench.gp comes from the fact that this program
is one (among many) examples where GP/PARI performs orders of magnitude
faster than systems such as Maple or Mathematica (try it!).

here's the script file:

u=1;v=1;p=1;q=1;
for(k=1,1000,w=u+v;u=v;v=w;p=p*w;q=lcm(q,w);if(k%10,,print(k,"
",log(p)/log(q))));


2) clareg.gp

Written entirely in the GP language without using the function buchgen,
the programs included in this file allow you in many cases to compute the
class number, the structure of the class group and a system of fundamental
units of a general number field (this programs sometimes fails to give an
answer). It can work only if initalg finds a power basis.
Evidently it is much less powerful and much slower than the function
buchgen, but it is given as an example of a sophisticated use of GP.
The first thing to do is to call the function clareg(pol,limp,lima,extra)
where pol is the monic irreducible polynomial defining the number field, limp
is the prime factor base limit (try values between 19 and 113), lima is
another search limit (try 50 or 100) and extra is the number of desired extra
relations (try 2 to 10). The program prints the number of relations that it
needs, and tries to find them. If you see that clearly it slows down too much
before succeeding, abort and try other values. If it succeeds, it will print
the class number, class group, regulator. These are tentative values. Then
use the function check(lim) (take lim=200 for example) to check if the value
is consistent with the value of the L-series (the value returned by check
should be close to 1). Finally, the function fu() (no parameters) returns a
family of units which generates the unit group (you must extract a system
of fundamental units yourself).

script:

rgcd(a,b,r)=a=abs(a);b=abs(b);while(b>0.01,r=a-b*trunc(a/b);a=b;b=r);a
{f(a,b,s,n,l)=s=a+b*t;n=abs(norm(s));nin=n;mv=vvector(li,j,0);
forprime(k=2,plim,l=0;while((n%k)==0,l=l+1;n=n/k);if(l,j=ind[k];cp=v[j][2];
while((a+b*cp)%k,j=j+1;cp=v[j][2]);mv[j]=l,));
if(n==1,lno=log(nin)/dep;vreg=vvector(lireg,j,if(j<=r1,log(abs(a+b*re[j])),log
(norm(a+b*re[j]))));
if(res,mreg=concat(mreg,vreg);m=concat(m,mv);
areg=concat(areg,a+b*t),mreg=mat(vreg);m=mat(mv);areg=[a+b*t]);
res=res+1;print1("(",res,": ",a,", ",b,")"),)}
{clareg(p,plim,lima,extra)=vi=initalg(p);p=vi[1];t=modp(x,p);
dp=disc(p);r=vi[6];r1=vi[2][1];findex=vi[4];
if(findex>1,print("sorry, the case f>1 is not implemented");1/0,);
print("discriminant = ",vi[3],", signature = ",vi[2]);
dep=length(p)-1;re=vector(lireg=(dep+r1)/2,j,if(j<=r1,real(r[j]),r[2*j-r1-
1]));
ind=vector(plim,j,0);v=[];
forprime(k=2,plim,w=lift(factmod(p,k));find=0;for(l=1,length(w[,1]),\
fa=w[l,1];if(length(fa)==2,if(find,,find=1;ind[k]=length(v)+1);
v=concat(v,[[k,-coeff(fa,0),w[l,2]]]),)));
li=length(v);co=li+extra;res=0;print("need ",co);
a=1;b=1;f(0,1);while(res<co,if(gcd(a,b)==1,f(a,b);f(-a,b),);
a=a+1;if(a*b>lima,b=b+1;a=1,));print(" ");mh=hermite(m);ms=matsize(mh);
if(ms[1]==ms[2],mhs=smith(mh);mh1=[mhs[1]];j=1;clh=mhs[1];
while(j<length(mhs),j=j+1;if(mhs[j]>1,mh1=concat(mh1,mhs[j]);clh=clh*mhs[j],))
;
print("class number = ",clh,", class group =
",mh1);km=kerint(m);mregh=mreg*km;
if(lireg==1,a1=1;print("regulator = 1"),coreg=length(mregh);
if(coreg<lireg-1,print("not enough relations for regulator: matrix size =
",matsize(mregh)),mreg1=extract(mregh~,cov=vector(lireg-
1,k,k))~;a1=det(extract(mreg1,cov));
for(j=lireg,coreg,a1=rgcd(a1,det(extract(mreg1,vector(lireg-1,k,k+j-
lireg+1)))));
print("regulator = ",a1))),print("not enough relations for class group: matrix
size = ",ms));}
{check(lim)=r1=vi[2][1];pol=vi[1];
z=2^(-r1)*(2*pi)^((r1+1-length(pol))/2)*sqrt(abs(vi[3]));
forprime(q=2,lim,z=z*(q-1)/q;fa=factmod(pol,q);
for(j=1,length(fa[,1]),z=z/(1-1/q^(length(fa[j,1])-1))));
clh*a1/rootsof1(vi)[1]/z}
{fu()=fw=[];for(k=1,length(km),ckm=km[,k];s=1;
for(j=1,length(ckm),s=s*areg[j]^ckm[j]);fw=concat(fw,s));fw}


3) lucas.gp

The half line function lucas(p) defined in this file performs the Lucas-Lehmer
primality test on the Mersenne number 2^p-1. If the result is 1, the Mersenne
number is prime, otherwise not.

4) rho.gp

A simple implementation of Pollard's rho method. The function rho(n) outputs
the complete factorization of n in the same format as factor.

script:

rho1(n,x,y)=x=2;y=5;while(gcd(y-
x,n)==1,x=(x*x+1)%n;y=(y*y+1)%n;y=(y*y+1)%n);gcd(n,y-x)
rho2(n,m)=m=rho1(n);if(isprime(m),print(m),rho2(m));if(isprime(n/m),print(n/m)
,rho2(n/m));
rho(n,m)=m=smallfact(n);print(m);n=m[length(m[,1]),1];if(isprime(n),,rho2(n));


5) squfof.gp

This defines a function squfof of a positive integer variable n, which may
allow you to factor the number n. SQUFOF is a very nice factoring
method invented in the 70's by D. Shanks for factoring integers, and is
reasonably fast for numbers having up to 15 or 16 digits. The squfof program
which is given is a very crude implementation. It also prints out some
intermediate information as it goes along. The final result is some factor of
the number to be factored.

script:

{squfof(n)=if(isprime(n),s=n,if(issquare(n),s=isqrt(n),p=smallfact(n)[1,1];
if(p!=n,s=p,if(n%4==1,dd=n;d=isqrt(dd);b=2*((d-1)\2)+1,dd=4*n;d=isqrt(dd);
b=2*(d\2));f=qfr(1,b,(b^2-dd)/4,0.);q=[];lq=0;ii=0;l=isqrt(d);flag=1;
while(flag,f=rhorealnod(f,d);ii=ii+1;a=compo(f,1);
if(ii%2,fl=0,fl=issquare(a));if(fl,as=isqrt(a);j=1;
while((j<=lq)&&fl,fl=(as!=q[j]);
j=j+1);if(as==1,s=0;flag=0,),);if(fl==0,if(abs(a)<=l,q=concat(q,abs(a));
print(q);lq=lq+1,),flag=0;gs=gcd(as,dd);print("i = ",ii);print(f);
gs=gcd(gs,bb=compo(f,2));
if(gs>1,s=gs,g=redrealnod(qfr(as,-bb,as*compo(f,3),0.),d);
fl=1;b=compo(g,2);while(fl,b1=b;g=rhorealnod(g,d);b=compo(g,2);fl=(b1!=b));
a=abs(compo(g,1));if(a%2,,a=a/2);s=a))))));s}


6) tutnf.gp

This is the sequence of GP instructions given in the tutorial in the section
on general number fields.

script:

\e
t=x^4+24*x^2+585*x+1791;nf=initalg(t);A=nf[1]
\precision=18
gc=galoisconj(A)
\precision=28
aut=gc[4]
pd=primedec(nf,7)
pr1=pd[1]
hp=idealmul(nf,idmat(4),pr1)
hp3=idealmul(nf,hp,idealmul(nf,hp,hp))
\\ or hp3=idealpow(nf,hp,3)
for(j=1,4,print(idealval(nf,hp3,pd[j])))
hpi3=[hp3,[0.,0.]];hr1=ideallllred(nf,hpi3,0)
hr=hr1;for(j=1,3,hr=ideallllred(nf,hr,[1,5]);print(hr))
arch=hr[2]-hr1[2];l1=arch[1];l2=arch[2];
s=real(l1+l2)/4;v1=[l1,l2,conj(l1),conj(l2)]~/2-[s,s,s,s]~
m1=nf[5][1];m=matrix(4,4,j,k,if(j<=2,m1[j,k],conj(m1[j-2,k])));
v=exp(v1);au=gauss(m,v)
vu=round(real(au))
u=mod(nf[7]*vu,A)
norm(u)
f1=factor(subst(char(u,x),x,x^2))
f1=factor(subst(char(-u,x),x,x^2))
v=sqrt(-v)
au=gauss(m,v)
v[1]=-v[1];v[3]=-v[3];au=gauss(m,v)
vu2=round(real(au))
u2=mod(nf[7]*vu2,A)
q=polred2(f1[1,1])
up2=modreverse(mod(q[3,1],f1[1,1]))
mod(subst(lift(up2),x,aut),A)
r=f1[1,1]%(x^2+u);-coeff(r,0)/coeff(r,1)
al=mod(x^2-9,A)
principalidele(nf,al)
for(j=1,4,print(j,": ",idealval(nf,al,pd[j])))
norm(al)
pd14=idealmul(nf,pd[1],pd[4])
idealmul(nf,al,idmat(4))
setrand(1);v=buchgenfu(A,0.2,0.2)[,1]
uf=mod(v[9][1],A)
uu=mod(v[8][2],A)
cu2=log(conjvec(u2));cuf=log(conjvec(uf));cuu=log(conjvec(uu));
lindep(real([cu2[1],cuf[1],cuu[1]]))
lindep([cu2[1],cuf[1],cuu[1],2*i*pi])
u2/uf
ru=nf[8]*vvector(4,j,coeff(v[8][2],j-1))
nf[7]*ru
setrand(1);bnf=buchinitfu(A,0.2,0.2);
bnf[8]
nf=bnf[7];
hp4=idealpow(nf,hp,4)
vis=isprincipal(bnf,hp4)
alpha=mod(nf[7]*vis[2],A)
norm(alpha)
idealmul(nf,idmat(4),mat(vis[2]))
vit=isprincipal(bnf,hp)
pp=isprincipal(bnf,pd14)
al2=mod(nf[7]*pp[2],A)
u3=al2/al
char(u3,x)
me=concat(bnf[3],[2,2]~)
cu3=principalidele(nf,u3)[2]
xc=gauss(real(me),real(cu3))
xd=cu3-me*xc
xu=principalidele(nf,uu)[2]
xd[1]/xu[1]
isunit(bnf,u3)
uu^3*uf



7) tutnfout

This is the slightly edited output of running tutnf.gp (obtained by removing
the ? at the beginning of each command for more legibility).

            GP/PARI CALCULATOR Version 1.38.23
                     (Sparcv8 version)

Copyright 1989,1993 by C. Batut, D. Bernardi, H. Cohen and M. Olivier

Type ? for help

\precision      = 28
\serieslength   = 16
\format         = g0.28
\prompt         = ?
stacksize = 4000000, prime limit = 500000, buffersize = 30000
? ? t=x^4+24*x^2+585*x+1791;nf=initalg(t);A=nf[1]
%1 = x^4 - x^3 - 21*x^2 + 17*x + 133
? \precision=18
   precision = 18 significant digits
? gc=galoisconj(A)
%2 = [x, 0, 0, -1/7*x^3 + 5/7*x^2 + 1/7*x - 7]
? \precision=28
   precision = 28 significant digits
? aut=gc[4]
%3 = -1/7*x^3 + 5/7*x^2 + 1/7*x - 7
? pd=primedec(nf,7)
%4 = [[7, [15, 12, 8, 8]~, 1, 1, [5, 4, 6, 4]~], [7, [7, 8, 8, 7]~, 1, 1, [5,
4, 2, 0]~], [7, [12, 13, 8, 7]~, 1, 1, [7, 2, 4, 0]~], [7, [13, 12, 8, 8]~, 1,
1, [5, 1, 4, 4]~]]
? pr1=pd[1]
%5 = [7, [15, 12, 8, 8]~, 1, 1, [5, 4, 6, 4]~]
? hp=idealmul(nf,idmat(4),pr1)
%6 =
[7 4 5 4]

[0 1 0 0]

[0 0 1 0]

[0 0 0 1]

? hp3=idealmul(nf,hp,idealmul(nf,hp,hp))
%7 =
[343 256 320 172]

[0 1 0 0]

[0 0 1 0]

[0 0 0 1]

? \\ or hp3=idealpow(nf,hp,3)
? for(j=1,4,print(idealval(nf,hp3,pd[j])))
3
0
0
0
? hpi3=[hp3,[0.,0.]];hr1=ideallllred(nf,hpi3,0)
%8 = [[7, 0, 6, 0; 0, 7, 2, 0; 0, 0, 1, 0; 0, 0, 0, 7], [-
4.542115783424510364590581431 - 4.894092298183431478504869242*i, -
3.241524812796742855830829541 - 6.227984987190759939849232263*i]]
? hr=hr1;for(j=1,3,hr=ideallllred(nf,hr,[1,5]);print(hr))
[[7, 0, 0, 3; 0, 7, 0, 1; 0, 0, 7, 5; 0, 0, 0, 1], [-
10.55920626489982283311310063 - 3.969550542422284695068789483*i, -
5.008074927542683607729721315 - 6.637335920436941617757515525*i]]
[[13, 0, 0, 10; 0, 13, 0, 2; 0, 0, 13, 11; 0, 0, 0, 1], [-
17.55370045383322168333975537 - 3.234391090803506441411369877*i, -
6.416260543236761408872612244 - 7.814752974898665176450527892*i]]
[[7, 0, 6, 0; 0, 7, 2, 0; 0, 0, 1, 0; 0, 0, 0, 7], [-
25.04390903221215492540062213 - 2.207432109741074676918714757*i, -
8.566810186297751680077560554 - 11.91890520832118510749679346*i]]
? arch=hr[2]-hr1[2];l1=arch[1];l2=arch[2];
? s=real(l1+l2)/4;v1=[l1,l2,conj(l1),conj(l2)]~/2-[s,s,s,s]~
%10 = [-3.794126968821658934140827421 + 1.343330094221178400793077242*i,
3.794126968821658934140827421 - 2.845460110565212583823780601*i, -
3.794126968821658934140827421 - 1.343330094221178400793077242*i,
3.794126968821658934140827421 + 2.845460110565212583823780601*i]~
? m1=nf[5][1];m=matrix(4,4,j,k,if(j<=2,m1[j,k],conj(m1[j-2,k])));
? v=exp(v1);au=gauss(m,v)
%12 = [79.00000000000000000000000006 + 9.693522799760103225000000000 E-27*i, -
24.00000000000000000000000002 + 1.137188718654215335000000000 E-27*i, -
14.00000000000000000000000000 - 2.909896133467555046000000000 E-28*i,
15.00000000000000000000000001 +  0.E-27*i]~
? vu=round(real(au))
%13 = [79, -24, -14, 15]~
? u=mod(nf[7]*vu,A)
%14 = mod(15/7*x^3 - 68/7*x^2 - 78/7*x + 79, x^4 - x^3 - 21*x^2 + 17*x + 133)
? norm(u)
%15 = 1
? f1=factor(subst(char(u,x),x,x^2))
%16 =
[x^8 + 85*x^6 + 1974*x^4 - 20*x^2 + 1 1]

? f1=factor(subst(char(-u,x),x,x^2))
%17 =
[x^4 + 13*x^3 + 42*x^2 - 8*x + 1 1]

[x^4 - 13*x^3 + 42*x^2 + 8*x + 1 1]

? v=sqrt(-v)
%18 = [0.09334880794728281557422432615 - 0.1174246256960511585313775710*i,
6.593348807947282815574224325 + 0.9834500294804898052951007426*i,
0.09334880794728281557422432615 + 0.1174246256960511585313775710*i,
6.593348807947282815574224325 - 0.9834500294804898052951007426*i]~
? au=gauss(m,v)
%19 = [-4.079490831526613599392561953 - 7.068193709477782249000000000 E-28*i,
0.9623261373956720055191059533 + 3.553714746609330177000000000 E-29*i,
1.007652680516380682499122747 + 9.093425419181585311000000000 E-29*i, -
0.9733355227802970462076159319 - 8.593538112938404083000000000 E-29*i]~
? v[1]=-v[1];v[3]=-v[3];au=gauss(m,v)
%20 = [-4.000000000000000000000000006 - 6.058451751247048377000000000 E-28*i,
1.000000000000000000000000002 - 8.884286869317293167000000000 E-29*i,
1.000000000000000000000000000 + 5.456055251881480216000000000 E-29*i, -
1.000000000000000000000000001 - 4.296769056469202041000000000 E-29*i]~
? vu2=round(real(au))
%21 = [-4, 1, 1, -1]~
? u2=mod(nf[7]*vu2,A)
%22 = mod(-1/7*x^3 + 5/7*x^2 + 1/7*x - 4, x^4 - x^3 - 21*x^2 + 17*x + 133)
? q=polred2(f1[1,1])
%23 =
[1 x - 1]

[1/7*x^3 + 2*x^2 + 7*x - 1/7 x^2 - x + 1]

[-x - 3 x^4 - x^3 - 21*x^2 + 17*x + 133]

[-4/7*x^3 - 7*x^2 - 21*x + 32/7 x^4 - 2*x^3 + 6*x^2 - 5*x + 133]

? up2=modreverse(mod(q[3,1],f1[1,1]))
%24 = mod(-x - 3, x^4 - x^3 - 21*x^2 + 17*x + 133)
? mod(subst(lift(up2),x,aut),A)
%25 = mod(1/7*x^3 - 5/7*x^2 - 1/7*x + 4, x^4 - x^3 - 21*x^2 + 17*x + 133)
? r=f1[1,1]%(x^2+u);-coeff(r,0)/coeff(r,1)
%26 = mod(1/7*x^3 - 5/7*x^2 - 1/7*x + 4, x^4 - x^3 - 21*x^2 + 17*x + 133)
? al=mod(x^2-9,A)
%27 = mod(x^2 - 9, x^4 - x^3 - 21*x^2 + 17*x + 133)
? principalidele(nf,al)
%28 = [[-9; 0; 1; 0], [4.071180332419999163417166456 -
2.351990513650678045072327170*i, -0.1793600343093725532064609704 +
1.836799691135712939544530673*i]~]
? for(j=1,4,print(j,": ",idealval(nf,al,pd[j])))
1: 1
2: 0
3: 0
4: 1
? norm(al)
%29 = 49
? pd14=idealmul(nf,pd[1],pd[4])
%30 =
[7 4 5 0]

[0 1 0 0]

[0 0 1 0]

[0 0 0 7]

? idealmul(nf,al,idmat(4))
%31 =
[7 4 5 0]

[0 1 0 0]

[0 0 1 0]

[0 0 0 7]

? setrand(1);v=buchgenfu(A,0.2,0.2)[,1]
%32 = [x^4 - x^3 - 21*x^2 + 17*x + 133, [0, 2], [18981, 7], [1, x, x^2,
1/7*x^3 + 2/7*x^2 + 6/7*x], [4, [4], [[7, 0, 0, 3; 0, 7, 0, 2; 0, 0, 7, 1; 0,
0, 0, 1]]], 3.794126968821658934140827422, 0.8826018286655581294913627961, [6,
1/7*x^3 - 5/7*x^2 - 8/7*x + 8], [1/7*x^3 - 5/7*x^2 - 1/7*x + 4], 150]
? uf=mod(v[9][1],A)
%33 = mod(1/7*x^3 - 5/7*x^2 - 1/7*x + 4, x^4 - x^3 - 21*x^2 + 17*x + 133)
? uu=mod(v[8][2],A)
%34 = mod(1/7*x^3 - 5/7*x^2 - 8/7*x + 8, x^4 - x^3 - 21*x^2 + 17*x + 133)
? cu2=log(conjvec(u2));cuf=log(conjvec(uf));cuu=log(conjvec(uu));
? lindep(real([cu2[1],cuf[1],cuu[1]]))
%36 = [0, 0, 1]
? lindep([cu2[1],cuf[1],cuu[1],2*i*pi])
%37 = [1, -1, -3, 0]
? u2/uf
%38 = mod(-1, x^4 - x^3 - 21*x^2 + 17*x + 133)
? ru=nf[8]*vvector(4,j,coeff(v[8][2],j-1))
%39 = [8, -2, -1, 1]~
? nf[7]*ru
%40 = 1/7*x^3 - 5/7*x^2 - 8/7*x + 8
? setrand(1);bnf=buchinitfu(A,0.2,0.2);
? bnf[8]
%42 = [[4, [4], [[7, 0, 0, 3; 0, 7, 0, 2; 0, 0, 7, 1; 0, 0, 0, 1]]],
3.794126968821658934140827422, 0.8826018286655581294913627961, [6, 1/7*x^3 -
5/7*x^2 - 8/7*x + 8], [1/7*x^3 - 5/7*x^2 - 1/7*x + 4], 150]
? nf=bnf[7];
? hp4=idealpow(nf,hp,4)
%44 =
[2401 942 1006 1544]

[0 1 0 0]

[0 0 1 0]

[0 0 0 1]

? vis=isprincipal(bnf,hp4)
%45 = [[0]~, [88, -31, -13, 15]~, 143]
? alpha=mod(nf[7]*vis[2],A)
%46 = mod(15/7*x^3 - 61/7*x^2 - 127/7*x + 88, x^4 - x^3 - 21*x^2 + 17*x + 133)
? norm(alpha)
%47 = 2401
? idealmul(nf,idmat(4),mat(vis[2]))
%48 =
[2401 942 1006 1544]

[0 1 0 0]

[0 0 1 0]

[0 0 0 1]

? vit=isprincipal(bnf,hp)
%49 = [[1]~, [43/7, -15/7, -6/7, 1]~, 139]
? pp=isprincipal(bnf,pd14)
%50 = [[0]~, [-40, 14, 6, -7]~, 143]
? al2=mod(nf[7]*pp[2],A)
%51 = mod(-x^3 + 4*x^2 + 8*x - 40, x^4 - x^3 - 21*x^2 + 17*x + 133)
? u3=al2/al
%52 = mod(-1/7*x^3 + 5/7*x^2 + 1/7*x - 4, x^4 - x^3 - 21*x^2 + 17*x + 133)
? char(u3,x)
%53 = x^4 - 13*x^3 + 42*x^2 + 8*x + 1
? me=concat(bnf[3],[2,2]~)
%54 =
[-3.794126968821658934140827422 + 10.76810805499055811618100739*i 2]

[3.794126968821658934140827422 + 6.579317850204167131564149548*i 2]

? cu3=principalidele(nf,u3)[2]
%55 = [-3.794126968821658934140827422 + 4.484922747810971639255720625*i,
3.794126968821658934140827422 + 0.2961325430245806546388627815*i]~
? xc=gauss(real(me),real(cu3))
%56 = [1.000000000000000000000000000,  0.E-48]~
? xd=cu3-me*xc
%57 = [ 0.E-47 - 6.283185307179586476925286766*i,  0.E-47 -
6.283185307179586476925286766*i]~
? xu=principalidele(nf,uu)[2]
%58 = [7.646841173435770929738625909 E-57 + 2.094395102393195492308428922*i, -
9.558551466794713662173282386 E-58 - 2.094395102393195492308428922*i]~
? xd[1]/xu[1]
%59 = -3.000000000000000000000000000 +  0.E-47*i
? isunit(bnf,u3)
%60 = [1, mod(3, 6)]
? uu^3*uf
%61 = mod(-1/7*x^3 + 5/7*x^2 + 1/7*x - 4, x^4 - x^3 - 21*x^2 + 17*x + 133)
?


----------------------- End example readme file --------------


I hope this helps to answer your question. :-)


--
Bye Manou
------------------------_----------------------------------------------
Manou BILLA       |  _ //  Connect your AMIGAs...
4, Ave. Nic Kreins|  \X/             ... A1000 / A2500 / A3000 ...
L-9536 WILTZ      |      ------ Member Team AMIGA Luxembourg ------
                  | email [e-mail]      FIDO 2:2455/560.8
-----------------------------------------------------------------------
           [PGP public key available on request]

... It is by coffee alone I set my mind in motion.
    It is by the beans of Java that thoughts aquire speed.
    The hands aquire shaking, the shaking is a warning.
    It is by coffee alone I set my mind in motion.